home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / rschur_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  1.6 KB  |  59 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    An upper quasi-triangular matrix.
  4.  
  5. // Syntax:    A = rschur ( N, MU, X, Y )
  6.  
  7. // Descriptioin:
  8.  
  9. //    A is an N-by-N matrix in real Schur form. All the diagonal
  10. //    blocks are 2-by-2 (except for the last one, if N is odd) and
  11. //    the k'th has the form [x(k) y(k); -y(k) x(k)]. Thus the
  12. //    eigenvalues of A are x(k) +/- i*y(k). MU (default 1) controls
  13. //    the departure from normality. 
  14.  
  15. //    Defaults: X(k) = -k^2/10, Y(k) = -k, i.e., the eigenvalues
  16. //                lie on the parabola x = -y^2/10.
  17.  
  18. //      References:
  19. //       F. Chatelin, Eigenvalues of Matrices, John Wiley, Chichester, 1993;
  20. //          Section 4.2.7.
  21. //       F. Chatelin and V. Fraysse, Qualitative computing: Elements
  22. //          of a theory for finite precision computation, Lecture notes,
  23. //          CERFACS, Toulouse, France and THOMSON-CSF, Orsay, France,
  24. //          June 1993.
  25.  
  26. //    This file is a translation of rschur.m from version 2.0 of
  27. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  28. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  29.  
  30. //-------------------------------------------------------------------//
  31.  
  32.  
  33. rschur = function ( n , mu , x , y )
  34. {
  35.   local ( n , mu , x , y )
  36.  
  37.   m = floor(n/2)+1;
  38.   alpha = 10; 
  39.   beta = 1;
  40.  
  41.   if (!exist (y)) { y = -(1:m)/beta; }
  42.   if (!exist (x)) { x = -(1:m).^2/alpha; }
  43.   if (!exist (mu)) { mu = 1; }
  44.  
  45.   A = diag( mu*ones(n-1,1), 1 );
  46.   for (i in 1:2*(m-1):2)
  47.   {
  48.     j = (i+1)/2;
  49.     A[i:i+1;i:i+1] = [x[j], y[j]; -y[j], x[j]];
  50.   }
  51.  
  52.   if (2*m != n)
  53.   {
  54.     A[n;n] = x[m];
  55.   }
  56.  
  57.   return A;
  58. };
  59.